Fit spatiotemporal tweedie models to biomass density and predict on grid
Max Lindmark
18 March, 2023
Intro
Fit Tweedie models to biomass density of cod, flounder, plaice and dab (juveniles and adults) between 1993-2020 in the Baltic Sea using sdmMTB fit different oxygen and temperature-related covariates. Predict on grid and use those predictions to calculate trends and velocities.
Load packages & source functions
# Load libraries, install if needed
library(tidyverse)
#> Warning: package 'tidyr' was built under R version 4.0.5
library(readxl)
library(tidylog)
library(RCurl)
library(devtools)
library(sdmTMB)
library(rMR)
library(kableExtra)
# Source code for plots
source_url("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/R/functions/map-plot.R")
#> Warning: package 'sf' was built under R version 4.0.5
theme_set(theme_plot())# Read data
d <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/data/clean/catch_clean.csv") %>%
rename(X = x, Y = y) %>%
drop_na(depth) %>%
mutate(depth_ct = depth - mean(depth),
depth_sc = depth_ct / sd(depth),
depth_sq = depth_sc^2,
year_f = as.factor(year),
quarter_f = as.factor(quarter)) %>%
pivot_longer(c(cod_adult, cod_juvenile, dab_adult, dab_juvenile, flounder_adult,
flounder_juvenile, plaice_adult, plaice_juvenile),
names_to = "group", values_to = "density") %>%
separate(group, into = c("species", "life_stage"), remove = FALSE)
#> Rows: 9376 Columns: 23
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (3): country, haul_id, ices_rect
#> dbl (20): year, lat, lon, quarter, month, sub_div, oxy, temp, sal, depth, x,...
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> rename: renamed 2 variables (X, Y)
#>
#> drop_na: removed 112 rows (1%), 9,264 rows remaining
#>
#> mutate: new variable 'depth_ct' (double) with 3,707 unique values and 0% NA
#>
#> new variable 'depth_sc' (double) with 3,707 unique values and 0% NA
#>
#> new variable 'depth_sq' (double) with 3,707 unique values and 0% NA
#>
#> new variable 'year_f' (factor) with 28 unique values and 0% NA
#>
#> new variable 'quarter_f' (factor) with 2 unique values and 0% NA
#>
#> pivot_longer: reorganized (cod_adult, cod_juvenile, dab_adult, dab_juvenile, flounder_adult, …) into (group, density) [was 9264x28, now 74112x22]
# Read metabolic parameter estimates and left_join
# mi_pars <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/data/clean/mi_params.csv") %>%
# dplyr::select(-...1, -temp)
# TODO: for now we'll use plaice parameters for flounder, see "00_estimate_mi_params.Rmd"
# mi_pars2 <- mi_pars
# mi_pars2[4, c(4:5)] <- mi_pars2[1, c(4:5)]
#
# d <- left_join(d, mi_pars2, by = "species")
# Read size csv to calculate the metabolic index
# sizes <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/data/clean/sizes.csv") %>%
# mutate(group = paste(species, name, sep = "_")) %>%
# dplyr::select(group, B)
#
# d <- left_join(d, sizes, by = "group")Read the prediction grid
# Read prediction grid
pred_grid <- bind_rows(readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/master/data/clean/pred_grid_(1_2).csv"),
readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/master/data/clean/pred_grid_(2_2).csv"))
#> Rows: 454412 Columns: 12
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (1): ices_rect
#> dbl (11): X, Y, year, lon, lat, depth, quarter, oxy, temp, sal, sub_div
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> Rows: 454412 Columns: 12
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (1): ices_rect
#> dbl (11): X, Y, year, lon, lat, depth, quarter, oxy, temp, sal, sub_div
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Calculate the metabolic index
# # Before we can calculate the pressure based metabolic index, we'll need to convert oxygen from ml/L mg/L and partial pressure
# # Some conversions
# # oxygen: https://www.ices.dk/data/tools/Pages/Unit-conversions.aspx
# # pressure: https://bluerobotics.com/learn/pressure-depth-calculator/
# # First calculate pressure in pascal: p_tot = (r*g*h) + p_atm, then to kgPa and finally atm
# rho = 1030 # density of water
# p_atm = 1.01325 # atmospheric pressure
# g = 9.807 # gravity of earth
#
# # Convert to partial pressure oxygen from ml/L
# d$po2 <- DO.unit.convert(d$oxy*(1/0.7), # d$oxy is in unit ml/L so we first need to convert to mg/L
# DO.units.in = "mg/L",
# DO.units.out = "PP",
# bar.press = ((rho*g*d$depth + p_atm)*10^-5)*0.986923,
# bar.units.in = "atm",
# bar.units.out = "kpa",
# temp.C = d$temp,
# salinity = 10, # TODO: d$sal, # let's fix salinity for easier comparison. It has a minor effect in this range
# salinity.units = "pp.thou")
#
# # Oxygen is ml/L, We want micro mol/L. 1 ml/l = 10^3/22.391 = 44.661 micro mol/l
# d$oxy_si <- (d$oxy * (10^3)) / 22.391
#
# # Calculate metabolic indices for a given mass and the temperature and oxygen in data
# # Line 123 in https://github.com/fate-spatialindicators/SDM_O2/blob/master/code/mi_functions.R
# kb <- 0.000086173324 # Boltzmann's constant
# t_ref <- 15 # arbitrary reference temperature
#
# # Calculate the metabolic index for pressure and concentration
# d <- d %>%
# mutate(inv_temp = (1/(temp + 273.15) - 1/(t_ref + 273.15)),
# phi = A0_po2*(B^n_po2)*po2 * exp((E_po2/kb)*inv_temp),
# phi_c = A0_o2*(B^n_o2)*oxy_si * exp((E_o2/kb)*inv_temp))Fit models and save AIC
# Fit models and save grid-level predictions
data_list <- list()
for(i in unique(d$group)) {
dd <- d %>%
filter(group == i) %>%
droplevels() %>%
drop_na(density)
mesh <- make_mesh(dd, c("X", "Y"), n_knots = 200)
plot(mesh)
m <- sdmTMB(
density ~ 0 + year_f + quarter_f + depth_sc + depth_sq,
data = dd, mesh = mesh, family = tweedie(link = "log"),
spatial = "on", time = "year", spatiotemporal = "IID",
control = sdmTMBcontrol(newton_loops = 1))
sanity(m)
# Plot residuals
res <- residuals(m) %>% as.data.frame()
ggplot(res, aes(sample = V1)) + stat_qq() + stat_qq_line() + theme(aspect.ratio = 1)
ggsave(paste0("figures/supp/qq_density_sdm/sdm_02_", i, ".pdf"),
width = 17, height = 17, units = "cm")
# Save model object
saveRDS(m, paste("output/models/sdm_02_", i, ".rds", sep = ""))
# Predict on grid
pred_grid <- pred_grid %>%
mutate(year_f = as.factor(year),
depth_ct = depth - mean(dd$depth),
depth_sc = depth_ct / sd(dd$depth),
depth_sq = depth_sc^2,
quarter_f = as.factor(quarter))
pred <- predict(m, newdata = pred_grid) %>% mutate(model = i)
data_list[[i]] <- pred
}
#> filter: removed 64,848 rows (88%), 9,264 rows remaining
#> drop_na: removed 13 rows (<1%), 9,251 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'year_f' (factor) with 28 unique values and 0% NA
#> new variable 'depth_ct' (double) with 11,144 unique values and 0% NA
#> new variable 'depth_sc' (double) with 11,144 unique values and 0% NA
#> new variable 'depth_sq' (double) with 11,144 unique values and 0% NA
#> new variable 'quarter_f' (factor) with 2 unique values and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 64,848 rows (88%), 9,264 rows remaining
#> drop_na: removed 13 rows (<1%), 9,251 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: changed 908,824 values (100%) of 'depth_ct' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sc' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sq' (0 new NA)
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 64,848 rows (88%), 9,264 rows remaining
#> drop_na: removed 629 rows (7%), 8,635 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: changed 908,824 values (100%) of 'depth_ct' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sc' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sq' (0 new NA)
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 64,848 rows (88%), 9,264 rows remaining
#> drop_na: removed 629 rows (7%), 8,635 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: changed 908,824 values (100%) of 'depth_ct' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sc' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sq' (0 new NA)
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 64,848 rows (88%), 9,264 rows remaining
#> drop_na: removed 235 rows (3%), 9,029 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: changed 908,824 values (100%) of 'depth_ct' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sc' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sq' (0 new NA)
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 64,848 rows (88%), 9,264 rows remaining
#> drop_na: removed 235 rows (3%), 9,029 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: changed 908,824 values (100%) of 'depth_ct' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sc' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sq' (0 new NA)
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 64,848 rows (88%), 9,264 rows remaining
#> drop_na: removed 622 rows (7%), 8,642 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: changed 908,824 values (100%) of 'depth_ct' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sc' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sq' (0 new NA)
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 64,848 rows (88%), 9,264 rows remaining
#> drop_na: removed 622 rows (7%), 8,642 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: changed 908,824 values (100%) of 'depth_ct' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sc' (0 new NA)
#> changed 908,824 values (100%) of 'depth_sq' (0 new NA)
#> mutate: new variable 'model' (character) with one unique value and 0% NA
# Save predictions and sims as data frames
df <- dplyr::bind_rows(data_list)
write_csv(df, "output/grid_pred_biomass_02.csv")